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We study collective plasmon excitations and screening of pure and disordered single- and bilayer 
black phosphorus beyond the low energy continuum approximation. The dynamical polarizability 
of phosphorene is computed using a tight-binding model that properly accounts for the band struc¬ 
ture in a wide energy range. Electron-electron interaction is considered within the Random Phase 
Approximation. Damping of the plasmon modes due to different kinds of disorder, such as resonant 
scatterers and long-range disorder potentials, is analyzed. We further show that an electric field ap¬ 
plied perpendicular to bilayer phosphorene can be used to tune the dispersion of the plasmon modes. 

For sufficiently large electric field, the bilayer BP enters in a topological phase with a characteristic 
plasmon spectrum, which is gaped in the armchair direction. 

PACS numbers: 73.21.-b, 73.22.Lp, 73.61.-r, 74.78.Fk 


I. INTRODUCTION 


Phosphorene is a new kind of two-dimensional (2D) 
material that can be obtained by mechanical exfolia¬ 
tion method from black phosphorus (BP) filmd^®. It 
is a semiconductor with direct band gap and highly 
anisotropic electronic and optical properties 3 The 

band gap of thin BP films varies, depending on the thick¬ 
ness, from 0.3 eV to around 2 eV. Furthermore, the elec¬ 
tronic band structure of this material is very sensitive to 
strain, 2 12 19 21 making it a promising candidate for elec¬ 
tromechanical applications. The plasmons and screen¬ 
ing in pristine single-layer and multilayer phosphorene 
were studied theoretically in Ref. [22] within a low en¬ 
ergy k • p model, where it was found that the dispersion 
of the collective modes in this material strongly depends 
on the direction of propagation. A similar approxima¬ 
tion was later used in Ref. (23l to study a generic double 
layer anisotropic system with large interlayer distance. 
Collective excitations of BP in the presence of a quan- 
tizying magnetic field have been studied recently, finding 
that the excitation spectrum is discretized into a series 
of anisotropic magneto-excitonsP^ 

A fundamental issue that needs to be addressed is the 
effect of disorder on the electronic and optical proper¬ 
ties of this material. Phosphorene and its la yered struc¬ 
tures are highly sensitive to the environment 5 due 

to its high reactivity when the samples are exposed to 
air. In this paper, we study the plasmon modes of dis¬ 
ordered single- and bilayer BP by using a tight-binding 
(TB) model which properly accounts for the electronic 
band dispersion in a wide energy window. 10 Damping 
of plasmons due to different kinds of disorder, such as 
point defects (resonant scatterers) and long-range disor¬ 
der potential, is analyzed. For this aim, the polariza¬ 
tion function is calculated with the tight-binding propa¬ 
gation method (TBPM^PMl 

which is extremely efficient 


in large-scale calculations of systems with more than mil¬ 
lions of atoms. The dielectric function is obtained within 
the random phase approximation (RPA). For the case of 
bilayer BP, we study the effect of an electric field applied 
perpendicular to the layers, which can be used to tune 
the dispersion of the plasmon modes. Finally we discuss 
the excitation spectrum of bilayer BP in the topologi¬ 
cal phase driven by the application of bias, in which the 
band structure is gapped in the armchair direction and 
dispersive in the zigzag direction. 

The paper is organized as follows. In Sec. [IT] we de¬ 
scribe the method used in our calculations. The dielec¬ 
tric screening properties are discussed in Sec. m and 
the plasmon excitation spectrum of disordered BP is an¬ 
alyzed in Sec. |IVl Our main conclusions are summarized 
in Sec. [VI 


II. MODEL AND METHOD 

BP is a single-elemental layered crystalline material 
consisting of only phosphorus atoms arranged in a puck¬ 
ered orthorhombic lattice. As sketched in Fig. [l] single 
layer phosphorene contains two atomic layers and two 
kinds of bonds with 2.22 A and 2.24 A bond lengths for 
in-plane and inter-plane P-P connections, respectively.® 
The electronic band structure around the gap can be de¬ 
scribed with a TB Hamiltonian for pristine BP with the 
forrrP^ 


U=T2 e i n i + Y + Y *jw c ! c i> (!) 

i i^j i^j 

where we consider five intralayer ( t\ = —1.220 eV, 1 2 = 
3.665 eV, t 3 = -0.205 eV, U = -0.105 eV, £5 = -0.055 
eV) and four interlayer hopping terms (t p 1 = 0.295 eV, 
t P 2 = 0.273 eV, t p s = —0.151 eV, t p 4 = —0.091 eV). 
For unbiased bilayer BP, there is an energy splitting of 
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Figure 1: Lattice structure of single-layer and bilayer BP. The 
relevant hopping terms considered in the Hamiltonian 0 are 
shown. The red and blue dots in (d) represent a resonant 
scatterer and a Gaussian center, respectively. 


Ae = 1.0 eV between the nonequivalent electrons in the 
sublayers. This model, which is based on first principle 
GW calculations, accurately reproduces the conduction 
and valence bands for energies ^0.3 eV beyond the gap. 
For a detailed description of the above TB Hamiltonian 
we refer the reader to Ref. urn 

We consider two main different sources of disorder 
in BP: local point defects and long range disorder po¬ 
tential. Point defects are modeled by elimination of 
atoms randomly over the whole sample, which can be 
viewed as vacancies, chemical adsorbates or substitution 
of other types of atoms which prevent the hopping of 
electrons to their neighbor sites P-^- 8 31 -^* This type of 
point defects are the so called resonant scatterers which 
provide resonances within the band gapl 18 * 35 * ^ , as con¬ 
firmed by fi rst-p rinciples calculations which include sin¬ 
gle vacancies 37 * 38 *, adatoms (Si, Ge, Au, Ti, V)^*, absorp¬ 
tion of organic moleculeJ^, substitutional p-dopants (Te, 
Cpor oxygen bridge-type defects 41 : 

The long-range disorder potential (LRDP), on the 
other hand, can account for electron-hole puddles, which 
are random changes of the local potential due to the inho¬ 
mogeneous distribution of charge in the sample. Within 
the TB model , it can be considered as a correlated Gaus¬ 
sian potential* 34 * 42 * 43 ^ such that the potential at site i fol¬ 
lows 

Vi = ^£4 exp j > ( 2 ) 

where r& is the position of the k -th Gaussian center, 
which are chosen to be randomly distributed over the 
centers of the projected lattice on the surface, £4 is the 
amplitude of the potential at the k-th Gaussian center, 
which is uniformly random in the range [—A, A], and d 
is interpreted as the effective potential radius. Typical 


values of A and d used in our model are A = 5 eV and 
d = 5a for LRDP.^ * 34 * 43 * Here a ~ 2.216A is the atomic 
distance of two nearest neighbors within the same plane. 
The origin of LRDP could be screened charged impu¬ 
rities on the substratd^® 7 or surface corrugations such 
as ripples and wrinkle^^ 9 . Although LRDP does not 
introduce resonances in the spectrum, they lead to an 
uniform increase of states in the gap™ The amount of 
resonant scatterers and Gaussian centers are quantified 
by n x and n c , respectively, which are the probability for 
a scatterer or Gaussian center to exist. 

The excitation spectrum and screening of single- and 
bilayer BP are calculated by using the TBPMpS® The 
dynamic polarization function is obtained from the Kubo 
formula^* 

j r°° 

n(q ,w) = —dre lu,T ([p(q,r),p(-q,0)]), (3) 

where V denotes the volume (or area in 2D) of the unit 
cell, p (q) is the density operator given by 


Nlayer 

p (q) = X! 52 c m c m ex P (*q • r ;,i) . ( 4 ) 

1=1 i 

and the average in ©> is taken over the canonical en¬ 
semble. For the case of the single-particle Hamiltonian, 
Eq. ([ 3 ]) can be written ad^ 

2 f°° 

H(q ,«;) = —-jf dr e iujT 

xlm (</j| n F ( H ) e lHr p (q) e~ lHr [1 - n F (ff)] p (-q) | p) , 

( 5 ) 

where rip {H) = ef3{H - M)+1 is the Fermi-Dirac distribu¬ 
tion operator, /3 = 1/kpT where kp is the Boltzmann 
constant and T is the temperature, and p is the chemical 
potential. We use units such that h = 1 and the average 
in Eq. © is performed over a random phase s nperp osi- 
tion of all the basis states in the real space, i.e.* 28 * 51 * 


\v) = 52 a M c Mi°)’ ( g ) 


where the coefficients dpi are random complex numbers 
normalized as J2u\ a i,i\ 2 = 1- We next introduce the 
time evolution of two wave functions 


ki(q^)> = e lHT [l-n F (H)]p(-(i)\ip ), (7) 

\<P 2 (r)) = e~ lHT n F {H)\ip), (8) 


which allows us to express the real and imaginary part 
of the dynamic polarization as 

2 f 00 

ReII(q,w) = -— dr cos(wr) Im (y> 2 ( T ) \p (q)| ¥>i ( T )) > 

2 f 00 

Imll (q, u>) = J dr sin(wr) Im (<£>2 (t)\p (q) \ { Pi ( r )) • 


( 9 ) 
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a matrix form 

\T( n \ — ( ^ ntra (^) Writer {q) A 

w V v intei (q) v intia (q) ) ’ 

where Vmtra(#) = V(o) as given by Eq. ( [TT] ), and the 
interlayer interaction between electrons in different BP 
layers is Einter = F(q)V(q) where F(q) = e~ qd and d = 
0.5455 nm is the interlayer distance. From this we can 
express the dielectric tensor for bilayer BP, including the 
intralayer and interlayer contributions, as 

£2 l (q,w) 

£il (q, uj) -V ( q ) F ( q ) n (q, uj) 

-V (q) F (q) n (q, uj) e 1L (q, u) 

( 12 ) 

In the next sections we use the polarization and dielec¬ 
tric functions defined here to study the effect of disorder 
on the screening properties and the collective plasmon 
modes of BP. To simplify the notation, we use e (q, uj) 
to represent £il (q, uj) for single-layer and the determi¬ 
nant of the dielectric tensor \£ 2 l (q, cj)| for bilayer. 



Figure 2: Static dielectric function of single- (left) and bilayer 
(right) BP. The chemical potential is set to be fi = 0.4 eV for 
single-layer and 0.73 eV for bilayer. The temperature is fixed 
as the room tempertuare T = 300 K. Plots (a) and (b) cor¬ 
respond to pristine BP, (c) and (d) correspond to samples 
with resonant scatterers originated from a concentration of 
n x — 0.2% vacant atoms, and (e) and (f) correspond to sam¬ 
ples with long range disorder potential with a concentration 
of n c = 0.2% Gaussian centers. 


The time evolution operator e~ lHr and Fermi-Dirac dis¬ 
tribution operator rip (H) can be obtained from the stan¬ 
dard Chebyshev polynomial decomposition. 28 ! In our cal¬ 
culations, the chemical potential is set to be ^ 0.1 eV 
above the edge of conduction band, i.e., /i = 0.4 eV for 
single-layer and 0.73 eV for bilayer BP. The temperature 
is fixed at T = 300 K. We use periodic boundary condi¬ 
tions, and the system size is 8192 x 8192 for single-layer 
BP, and 6000 x 6000 for bilayer. 

Electron-electron interactions are considered within 
the RPA. The two-dimensional dielectric function for 
single-layer phosphorene is calculated as 

SiL (q, w) = 1 - V (q) n (q, w) (10) 

where 

<“> 

is the Fourier component of the Coulomb interaction in 
two dimensions, in terms of the background dielectric 
constant k. For bilayer BP, the Coulomb interaction has 


III. DIELECTRIC SCREENING 

We start by studying the effect of disorder on the 
screening properties. In Fig. [2j we plot the static di¬ 
electric function 5 (q, 0) of single- and bilayer BP. There 
is a clear dependence of the static dielectric function with 
the direction of momentum, showing a different behavior 
for q along the zigzag (X) or armchair (Y) directions. 
For the pristine case, our numerical results for the mo¬ 
mentum dependence of the static dielectric function of 
single-layer phosphorene, Fig. [2|a), agree with the ana¬ 
lytic low energy model presented in Ref. [221 The present 
calculation within the TBPM method allows us to extend 
the calculation of Ref. 22| to larger wave-vectors and to 
consider realistic impurities in the sample. We notice 
that the size of the lattices used in our simulations deter¬ 
mines the rage of wave-vectors for our numerical results, 
which is q > 0.01a _1 . We do not perform calculations 
with momentum transfer smaller than q = 0.01a -1 since, 
in order to preserve the numerical accuracy for q —)> 0, 
an extreamly large sample size N oo is required. 

In both single- and bilayer BP, point-like resonant scat¬ 
terers lead to the creation of midgap states inside the 
gap. 18 These states are quasi-localized around the cen¬ 
ter of the impurity. Therefore, the typical divergence of 
the dielectric function at q 0, characteristic of metallic 
behavior, is reduced for BP crystals with this kind of dis¬ 
order, as it can be seen in Fig. 2jc)-(d). A similar result 
for the dielectric function is obtained for LRDP which 
can be associated to the existence of electron-hole pud¬ 
dles in the sample. A drop of the dielectric function at 
small q appears also in disordered graphene with chem¬ 
ical absorbers. 30 The effect of disorder on the static di¬ 
electric function is negligible at short wave-lengths (large 
q vectors). 
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Figure 3: Energy loss functions of single-layer BP along zigzag 
(left) and armchair (right) directions. The simulations are 
done for: (a,b) pristine BP, (c,d) samples with resonant scat¬ 
tered, and (e,f) samples with long-range disorder potentials. 
All the results are calculated at T — 300 K. We notice that the 
apparent discretization of the spectrum in some of the plots 
is an artifact due to finite size limitations in our calculations. 


IV. ENERGY LOSS FUNCTION AND 
PLASMONS 

In this section we study the collective modes of the 
system, which are defined by the zeroes of the dielec¬ 
tric function £ (q, uS) for single-layer BP (or zeroes of the 
dielectric tensor determinant for bilayer BP). The dis¬ 
persion relation of the plasmon modes is defined from 
Re e(q ,u p i) = 0, which leads to poles in the energy loss 
function — $sl/e(q, uj), that can be measured by means of 
electron energy loss spectroscopy (EELS)f^ 


A. Single-layer BP 

For single-layer BP, our numerical results for the loss 
function are plotted in Fig. [3] and show anisotropic plas¬ 
mon modes along the zigzag and armchair directions, 
both following a low energy y/q dependence at small mo¬ 
mentum transfer. This is due to the paraboloidal band 
dispersion of BP at low energies, and it is consistent with 
the results for the excitation spectrum of BP within the 
k • p approximation. 22 The present method allows us to 
study the spectrum at higher energies and wave-vectors, 
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Figure 4: Energy loss functions of single-layer BP along zigzag 
(left) and armchair (right) directions. Panels (a) and (b) com¬ 
pare the loss function of pristine single-layer BP with that of 
samples with different concentrations of point defects. Panels 
(c) and (d) correspond to LRDP. All the results are calculated 
at T = 300 K. 


including also the effect of disorder in the sample. 

In the presence of disorder, there is a broadening and 
red shift of the plasmon modes, as it can be seen in Fig. 
|3jc)-(d) for single-layer samples with resonant scatterers, 
and in Fig. |3je)-(f) for long range disorder potential. For 
a fixed wave-vector, the effect of disorder is better seen 
by the loss function shown in Fig. [4] for q = 0.1a _1 
and different concentrations of disorder. Notice that the 
presence of LRPD induces a stronger damping of the 
plasmon, as compared to point defects. This is seen by 
comparing the broadening of the peaks of the top and 
bottom panels of Fig. [4j Since plasmons are collective 
electronic oscillations induced by long range Coulomb in¬ 
teraction, it is expected that long range disorder affects 
more the coherence of these modes than local scatterers 
as vacancies, which are relevant at short length scales. 

Fig. @ also shows that, for a given wavelength, the 
two kinds of disorder considered in our simulations lead 
to a shift of the plasmon resonance towards lower en¬ 
ergies. This effect can be understood from a perturba¬ 
tive point of view, by using the disordered averaged re¬ 
sponse function introduced by Mermin in the context of a 
2DEG with a parabolic band dispersion,^ and which has 
been recently used to study the effect of weak disorder on 
the plasmon dispersion of spin-polarized grapheneP^ In 
the weak disorder limit, and assuming for simplicity an 
isotropic parabolic band dispersion, the disordered aver¬ 
aged response function has the form 53 55 




~N( 0) 1 


“ _ ) 

V(w + i/TeO 2 -v|<7 2 - */ r eZ / 


( 13 ) 

where r e i is the elastic life-time of the momentum eigen¬ 
states of the disordered system, N(0) = grrib/2'K is the 
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density of states (per unit area) at the Fermi energy, 
where g = 2 is the spin degeneracy, is the effective 
mass, and is the Fermi velocity. In the presence of 
disorder the dispersion relation of the plasmon will have, 
even in the low energy and long wavelength limit, and 
imaginary part, due to finite damping of the plasmon. 
Within this approximation, it is possible to obtain an 
analytical expression for the plasmon dispersion by sub¬ 
stituting (13) into ( 10 ) and then finding the zeros of the 
dielectric function, with the approximate solution 


UpM) « g + (y^(<? + 8r sM - 1/7 \ - i/Tei'j 

(14) 

where r s = 2 m^e 1 / nkp is the dimensionless interaction 
parameter of a 2DEG and kp is the Fermi wave-vector. 
Eq. ( p~4| ) reduces to the standard plasmon dispersion in 
the clean (1 /r e i 0) limit. We notice here that the 

presence of impurities limit the dispersion of collective 
plasmon modes for length scales longer than the mean 
free path, defining a minimal wave-vector below which 
the plasmon mode is not well defined. The existence 
of such length scale has been discussed, in the context 
of graphene, in Ref. |5U For the case of interest here 
su ch infrared cuto ff takes the form g min = —Ar s kp + 
V 78 rjk% + l/T&vjr, below which plasmon modes do not 
exist. 

The next step in our analysis is to quantify the effect 
of disorder on plasmon losses. The physical observability 
of a damped collective mode depends on the sharpness of 
the resonance peak, and it is defined as a dimensionless 
inverse quality factor 7Z = 7 /uj p i, where the damping rate 
7 of the mode, which is proportional to Im n(q,cOp/), is 
given by 


7 = 


Im II(q, Wpi) 
£ReII(q,w)7 


(15) 


The plasmon lifetime r is calculated via the inverse 
quality factor as r = \/'Rjuj p i. We have calculated nu¬ 
merically 7 Z ( 7 ) and r, and the results are shown in Fig. 
[5] We first notice that, for a given concentration of im¬ 
purities and for a given wavelength, the plasmon rate is 
anisotropic. In fact, our results show that plasmons dis¬ 
persing in the ^/-direction (armchair) are more efficiently 
damped than plasmons dispersing in the x (zigzag) direc¬ 
tion, which have a larger lifetime. Second, as discussed 
above, we find that short range resonant scatterers like 
point defects cause less losses in the plasmon coherence 
than LRDP. 


B. Bilayer BP 

In the following we present results for bilayer BP, con¬ 
sidered as two phosphorene layers coupled by interlayer 
hopping (see Fig. [I]), and interacting via long range 
Coulomb potential, as explained in Sec. [IlJ Our results 
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Figure 5: Effect of disorder on plasmon losses. Panels (a) and 
(b) show the plasmon life-time, and panels (c) and (d) show 
the corresponding damping rates (inverse quality factor). The 
results for bilayer BP samples correspond to the in-phase ~ 
yfq plasmon mode (see text). 


for the loss function are shown in Fig. [ 6 ] The first thing 
to notice is the existence of two collective modes in the 
spectrum, as it is usual in bilayer systems. 23 * 56 ! ^ One 
is the classical cj+(g —» 0 ) ^ y/q plasmon mode, which 
has its counterpart in single-layer samples (see Fig. [ 3 ]). 
This mode correspond to a collective excitation of the 
electron liquid in which the carriers of both layers os¬ 
cillate in-phase. In addition to this plasmon, in Fig. [ 6 ] 
we observe the existence of an additional mode in the 
excitation spectrum dispersing below the plasmon. 
This mode, with a low energy linear dispersion relation, 
(q 0) ~ g, is characteristic of bilayer 2DEG systems 
and corresponds to a collective oscillation in which the 
carrier density in the two layers oscillates out-of-phase. 
The existence of this mode is also clear in Fig. [7| where it 
can be seen that, besides the peak in the loss function cor¬ 
responding to the in-phase plasmon, there is a second 
resonance at lower energies (with less spectral weight) 
which corresponds to the ^ q plasmon (indicated by 
black arrows). The uj- mode is much weaker as compared 
to the mode because of the short distance between the 
two layers, which complicates the out-of-phase density os¬ 
cillation. In Ref. [23j the two modes discussed above have 
been also obtained by using a general low energy model 
of anisotropic double-layer systems. We notice here that, 
since they study bilayer systems in which the two layers 
are separated by more than 5 nmp 3 their results show a 
lO- mode that is more coherent and dispersing than the 
one obtained in our calculations. As a matter of fact, 
the dispersion relation of the out-of-phase mode is pro¬ 
portional to the layer separation, cj_(g) ^ d}' 2 q , where 
d is the interlayer distance. Here we restrict our study 
to real bilayer phosphorene, in which the two layers are 
separated by d = 0.524 nm, and interlayer hopping of 
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Figure 7: Energy loss functions of bilayer BP along zigzag 
(left) and armchair (right) directions. Panels (a) and (b) com¬ 
pare the loss function of pristine bilayer BP with that of sam¬ 
ples with different concentrations of point defects. Panels (c) 
and (d) correspond to LRDP. All the results are calculated 
at T — 300 K. The peak marked by arrows in each panel 
correspond to the out-of-phase cj~ ( q ) ~ q plasmon (see text). 


Figure 6: Energy loss functions of bilayer BP along zigzag 
(left) and armchair (right) directions. The simulations are 
done for: (a,b) pristine bilayer BP, (c,d) samples with res¬ 
onant scatterers, and (e,f) samples with long-range disorder 
potentials. All the results are calculated at T — 300 K. 


electrons is allowed. As a consequence, the uj- mode ob¬ 
tained in our simulations for bilayer BP disperses with a 
smaller velocity than the mode of Ref. [23]. 

We also observe that, as in the single-layer case, both 
modes are damped due to disorder. The characteris¬ 
tics of the damping of the in-phase uj+ ~ ^Jq mode can 
be understood from the previous analysis for single-layer 
BP. Losses are stronger in the out-of-phase uj- plasmon, 
whose coherence basically vanishes due to disorder. The 
reasons for this are twofold. On the one hand, the spec¬ 
tral weight of the mode is considerably weaker than 
the plasmon. On the other hand, the uj- mode dis¬ 
perses in the uj — q spectrum very close to the electron- 
hole continuum, whose threshold can be reached due to 
disorder (and thermal) broadening of the plasmon peak, 
with the consequent Landau damping of the mode. 


C. Biased Bilayer BP 

In this section we include in our calculations the pres¬ 
ence of a perpendicular electric field applied to the bi¬ 
layer BP. This technique can be used to manipulate the 
electronic band structure of the system, and it has been 
proposed recently as an appropriate way to drive a nor¬ 
mal to topological phase in this material. 59 60 Here we 


are interested on studying the effect of an electric field 
on the dispersion relation of plasmons. For this aim, we 
introduce a biased on-site potential difference A between 
the two layers as described in Ref. 61, and consider three 
representative cases, as shown in Fig. |8j A = 1 eV, for 
which the band gap is reduced, A = 1.4 eV, for which 
the gap completely closes, and A = 1.8 eV, for which 
the conduction and valence bands are overlapped, cor¬ 
responding to the topological phase discussed in Refs. 
1591601 It is interesting to notice that, for this last case, 
the band structure of bilayer BP presents Dirac like cones 
for the dispersion in the k x direction [Fig. |8](e)jL and the 
spectrum is gapped in the k y direction [Fig. p[f)]. For 
A = 1 eV and 1.4 eV, the chemical potential is chosen to 
be about 0.1 eV above the edge of the conduction band, 
to better compare with the results of unbiased bilayer BP 
of Fig. [6j 

Our results show that the velocity of the plasmon mode 
increases with the biased voltage, as it can be seen from 
the evolution of the excitation spectrum with the applied 
A in Fig. [9j This suggest the possibility to tune the 
plasmonic properties of this material with a perpendicu¬ 
lar electric field. Interestingly, when the applied electric 
field exceeds the critical field and enters in the topologi¬ 
cal phase, the plasmons are more coherent, as it can be 
inferred by looking at the strength of the modes in the 
density plots of Fig. [9] (notice the different scales in the 
density bar of each panel). Furthermore, the plasmon 
dispersion is ungapped along q x direction, whereas the 
presence of a gap in the electronic band structure in the 
T—Y direction leads to a gapped collective mode along 
q y , as it can be seen by looking Fig. [9](f ) in the low en¬ 
ergy and small wave-vector region. We notice again that, 
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Figure 8: Low energy band structure of bilayer BP in the 
presence of a perpendicular electric field, for different biased 
potential energy between top and bottom layers. The red 
dashed line in each panel indicates the position of the chemical 
potential in the calculation of the plasmon modes. 

due to the finite size nature of our simulations, our re¬ 
sults cannot reach the q 0 limit, for which we would 
need to consider an infinite sample. However, the exis¬ 
tence of a gap in the plasmon spectrum q y direction is 
already clear in the present calculation for A = 1.8 eV. 
Notice that the peculiar band structure in the topologi¬ 
cal phase, with Dirac like band crossings in one direction 
and gap in the other direction, leads to a rich excita¬ 
tion spectrum with a numerous peaks corresponding to 
the plasmon collective excitations and to the enhanced 
optical transitions due to appearance of Van Hove singu¬ 
larities in the spectrum. 


V. CONCLUSIONS 

In conclusion, we have studied the effect of disorder 
in the excitation spectrum of single-layer and bilayer BP. 
The band structure has been calculated with an accurate 
tight-binding model which is valid within ~ 0.3 eV be¬ 
yond the gap. The dynamical polarization function has 
been calculated with the Kubo formula, and from this, 
the energy loss function has been obtained within the 
RPA. We have found that disorder leads to a redshift of 
the plasmon resonance. This effect has been discussed 
from a perturbative point of view, within the framework 
of the disordered averaged response function. 53 
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Figure 9: Energy loss functions of biased bilayer BP along 
zigzag (left) and armchair (right) directions. The strength of 
the biased potential is given in each panel. 


The different kinds of disordered that have been ana¬ 
lyzed show that, for the same concentration of impuri¬ 
ties, LRPD leads to stronger damping rates as compared 
to local point defects. Such effect is understood from 
the intrinsic long wavelength nature of the plasmon os¬ 
cillations, induced by long range Coulomb interaction. 
Therefore resonant scatterers as point vacancies, which 
modify the electronic properties at very short length 
scales, are less effective inducing losses of the plasmon 
modes than LRPD. 

The spectrum of bilayer BP presents two plasmon 
modes. One mode in which the carrier density in the two 
planes oscillates in-phase, with a dispersion cj + (g) ~ yfq, 
which is the counterpart of the standard plasmon in 
single-layer BP, and one mode in which the carriers in 
the two layers oscillate out-of-phase, with a liner disper¬ 
sion relation at low energies, c o~(q) ~ q . The coherence 
of the cj_ mode is more dramatically affected by disorder. 

Finally, we have studied the effect of a perpendicu¬ 
lar electric field in the excitation spectrum of bilayer 













































BP. We have shown that the dispersion of the collective 
modes can be tuned by the application of such biased 
field. Furthermore, we have shown that, beyond some 
critical field, the bilayer BP enters in a topological phase 
with Dirac like crossing of bands in one direction and 
gapped in the other direction. As a consequence, the ex¬ 
citation spectrum is highly rich in this limit, with highly 
coherent plasmon modes which are gapped in the q y di¬ 
rection of the spectrum. In summary, our results show 
a highly anisotropic excitation spectrum of single layer 
and bilayer BP, features that could be of high interest 
for future optoelectronic applications. 
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